Self-interaction errors in density functional calculations of electronic transport 
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All density functional calculations of single-molecule transport to date have used continuous 
exchange-correlation approximations. The lack of derivative discontinuity in such calculations leads 
to the erroneous prediction of metallic transport for insulating molecules. A simple and com- 
putationally undemanding atomic self-interaction correction greatly improves the agreement with 
experiment for the prototype Au/dithiolated-benzene/Au junction. 
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Molecular devices are becoming increasingly important 
in a wide spectrum of applications. These range from the 
building blocks of revolutionary computer architectures 
0, to disposable electronics, to diagnostic tools for ge- 
netically driven medicine, to multifunctional sensors |2j. 

Therefore interest is growing in the development 
of computational tools capable of predicting the I-V 
characteristics of devices comprising only a handful of 
atoms. In general these are based on Landauer scat- 
tering theory, typically in the non-equilibrium Green's 
function (NEGF) formalism combined with an elec- 
tronic structure method, usually density functional the- 
ory (DFT) 0, . Such schemes are physically appealing 
andyield useful results Q, even if they are incomplete 
0, 0- an d are computationally simpler than many- 
body methods [Tof . The fundamental requirements for 
an electronic structure theory applied to the problem of 
transport through single molecules are: 1) to be accurate 
when the molecule has a fractional number of electrons, 
2) to describe properly both the electron affinity (A) and 
ionization potential (I) of the isolated molecule, 3) to be 
cast in a single particle form. The first two conditions 
are necessary for a correct description of the transport, 
while the third produces computationally efficient algo- 
rithms. Most importantly, as we show here, we need an 
accurate description of the HOMO state as a function of 
its occupation. 

The exact Kohn-Sham (KS) potential of a ^-electron 
system always satisfies the condition Chqmo = ~~In, i-e., 
the highest occupied KS molecular orbital energy is the 
negative of the ^-electron ionization potential. Let N+n 
be the number of electrons localized on a molecule weakly 
coupled to a reservoir, where N is an integer, but n is 



< n < 1, e 
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To achieve this, the KS 



potential jumps by a step of In — ijv+i = In — An, 
where An is the electron affinity. This is the infamous 
derivative discontinuity of DFT 0, j which is miss- 
ing in ordinary continuous functionals such as the local 
density (LDA) or the generalized gradient approxima- 
tion (GGA) |l3|. Smooth exchange-correlation function- 
als continuously connect the orbital levels for different 
integer occupations, leading to qualitative errors such as 
the erroneous prediction of the dissociation of heteronu- 
clear molecules into fractionally charged ions We 
now show the errors that this gives rise to in a typical 
transport calculation. 

We model a two terminal molecular device as two 
featureless leads (constant density of states) kept at 
different chemical potentials /il and /ir and coupled 
through a single energy level e (Fig. IHa)). The den- 
sity of states (DOS) associated with e is a Lorentzian, 
^C^O = 2tt (E-e)*+(r JW ' w ith broadening T arising from 
the hopping to the leads. The energy level occupation n 
and the steady state current / can be obtained by bal- 
ancing the in-going and out-going currents to and from 
the energy level |3j. At steady state, n is just propor- 
tional to the Fermi distributions f(e,T) of the leads: 
n = dE D(E)[f(E - /i L , T) + f(E - m , T)], while 
the current is given by 



7 =^ 



dED(E)[f(E- f i L ,T)-f(E- m ,,T)} . (1) 



continuous. For — I < n < 0, e 
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The dynamics of this model are rather simple. Con- 
sider the weak coupling limit (r An), where simply 
D(E) ~ S(E—e). Both occupation and current are solely 
determined by the position of the energy level with re- 
spect to the chemical potential of the leads. If e is larger 
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than both fii, and /iR, then n « and no current flows. 
In contrast, if the energy level is below the chemical po- 
tentials of both leads, then n w 2, but the current is still 
zero. Finally, if /ir < e < Ml the occupation is < rt < 2 
and current flows. Considering now that the chemical 
potential in the leads is simply Ml/r = E f ± eV/2, where 
E F is the Fermi level of both leads and V is the applied 
bias, this simple model predicts a conductance gap in the 
I-V curve for -2\E F - e| < eV < 2\E F - e\. 
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FIG. 1: (a) Schematic two terminal device. Two leads are 
kept at the chemical potentials /iL and /ir and the transport 
is through a single energy level e. The hopping energy between 
the leads and the energy level is V. (b) Dependence of e on its 
own occupation n. The straight line corresponds to a typical 
LDA dependence and the step-like line to the DISC. Notice 
that the LDA line becomes flat at n ~1.5 (e(n)=0) since here 
the eigenstate is unbound. The dotted horizontal line denotes 
the position of the leads Fermi level Ef. 

However, because this is an effective one-body repre- 
sentation of an interacting system, in general the posi- 
tion of the energy level depends on its own occupation, 
e — e(n). Let us now solve this problem within KS DFT 
For defmiteness, assume that e is the LUMO (n=0) 
of a certain molecule, which contains N electrons in the 
neutral state. In the exact KS theory, when this molecule 
is weakly coupled to a reservoir, e will be a discontinu- 
ous function of n 14] . We parameterized this depen- 
dence with the step-like curve of figure ^b), and we call 
it "DISC" (discontinuous occupation). For < n < 1, 
e(n) = — An, where An is the electron affinity of the 
isolated molecule, while e(n) jumps rapidly to its next 
plateau (— -Ajv+i) just above 1. In contrast, the LDA en- 
ergy level position e varies approximately linearly with 
n (e = U n), reflecting the fact that the LDA total en- 
ergy varies approximately quadratically around the neu- 
tral configuration (lj) . 

The different I-V characteristics that one obtains by 
using either the LDA or DISC parameterizations are pre- 
sented in Fig. [3 along with the level occupation and its 
position as a function of bias. These have been obtained 
by iterating self-consistently n(e) with e(n), where we 
assume e(0) just below E F (|e(0) — E F \ = 0.5 cV). Con- 
sider first the weak coupling limit. In both LDA and 
DISC, the energy level pins E F at zero bias. As the bias 
is further increased, more charge fills the energy level, 
which keeps rising up. Figure |2Icl) shows that this rise 



is found both in LDA and DISC and is approximately 
linear with the bias. Importantly, as soon as e shifts 
above the chemical potential of the right-hand side con- 
tact, then /(e — /l«r) 0, and the current will be simply 
proportional to the level occupation (J « T/(e — /Xl)). 

Clearly LDA and DISC behave in a qualitatively differ- 
ent way. In fact, a LDA- type potential leads to a linear 
dependence of the occupation on bias (see Fig. |2Jbl)), 
and consequently to a metallic conductance. In contrast, 
in DISC, the energy level shifts upwards without substan- 
tial charging. The result is that the occupation jumps 
almost discontinusly from n — to n — 1 when the bias 
is increased, and consequently a gap in the I-V curve 
opens. Such a gap is as large as the one in the occupa- 
tion, which is roughly ^4^. 

We emphasize that, despite the simplicity of the func- 
tion used to mimic the discontinuity, our model restores 
the correct I-V behavior with the expected conductance 
gap, roughly equal to the A jy, thus repairing the faulty 
LDA description. The lack of eigenvalue discontinuity 
causes a dramatic overestimation of the metallicity for 
a molecular junction obtained within LDA. This result 
generalizes arguments previously given using many-body 
theory 0, or, within DFT, only for weak bias 0. 
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FIG. 2: (a) Current I, (b) occupation n and (c) position of the 
energy level e as a function of bias V. The parameters used 
here are e(0)=-5.5 eV, (7=5 eV, E F =-5.0 eV and T=300°K. 
The curves on the left-hand side are obtained in the weak 
coupling limit (r=0.2 eV) and those on the right-hand side 
in the strong coupling limit (r=1.2 eV). 

In contrast, the curves obtained for LDA and DISC in 
the strong coupling limit (figures [3^a,2) , EJb2) and[2Ic2)) 
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look rather similar. This is because a considerable frac- 
tion of the level occupation and the current comes from 
the tail of the energy level DOS. In the large coupling 
limit, r ~ Ajv, both n and the current I are rather in- 
sensitive to e(n), and we find that standard continuous 
functionals give rather accurate I-V characteristics. 

Having identified the lack of the derivative disconti- 
nuity in LDA as a major source of error in DFT-based 
transport calculations, we propose a corrective scheme 
for the NEGF method. The key consideration is that 
in LDA (or GGA), the linear behavior of the KS eigen- 
values and the absence of derivative discontinuity in the 
total energy functional, is mainly due to the presence of 
self-interaction error (SIE), that is, the interaction of an 
electron with the exchange and correlation potential gen- 
erated by its own charge This spurious interaction 
is responsible for a series of failures of DFT. Most no- 
tably, negatively charged ions are unstable and in general 
~ £ homo ^ s n °t even close to the ionization potential. The 
elimination of the SIE improves considerably the agree- 
ment with experiments with respect to LDA and, more 
importantly in this context, makes the KS eigenvalues 
resemble more closely the true removal energies [l6| . 

Direct subtraction of SI in atoms is conceptually and 
practically simple However, the application of the 

method to extended systems is both cumbersome and 
computational demanding A useful alternative is 

that first proposed by Vogel and then extended by Fil- 
ippetti in which the SI is parametrized in terms of its 
atomic counterparts and subtracted out (Pseudo-SIC, 
PSIC) [lfjj . In the spirit of this method we have con- 
structed an effective tight-binding model, and investi- 
gated the transport of a benzene-(l,4)-dithiolate (BDT) 
molecule sandwiched between two (001) oriented gold 
current/ voltage probes (Fig.[3J). This represents a bench- 
mark, since LDA-DFT calculations an d experiments 
[2flj are in stark disagreement. In particular most of the 
calculations fail to predict a conductance gap at zero bias. 

Here we adopt a minimal ir model where we consider 
only p z orbitals (orthogonal to the BDT plane) for both 
C and S atoms and s orbitals for Au. H atoms are sim- 
ply used for passivation and are not considered explic- 
itly. The on-site energy of such p z orbitals are parame- 
terized from their atomic counterparts and coincide with 
the HOMO state of the free atom. This is computed for 
different occupations with a standard self-consistent cal- 
culation using either LDA and SIC 0|. The resulting 
e(n) curves (not reported here) are similar to that of fig- 
ure I2b). Our procedure neglects the crystal field, and 
assumes that the electron screening is weak. Although for 
a fully quantitative analysis such aspects must be consid- 
ered, we do not expect that these details will change the 
main features of our model. Finally, the hopping inte- 
grals are taken from the literature [2 1| . 

The I-V characteristics are then calculated using stan- 
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FIG. 3: BDT molecule attached to (001) oriented Au surfaces. 
The angle cf> between the BDT plane and the gold modulates 
the strength of the molecule/lead coupling. 



dard NEGF method olog y with a tight-binding version of 
our code Smeagol [6ll22|. In the simulation, we alter the 
strength of the coupling to the leads by varying the angle 
(j) between the BDT plane and the apex of the Au pyra- 
mid (see figure|3J) . The coupling is then 7 sin <fr with 7 the 
Au-S spa hopping integral. The alignment of E-p of the 
leads with chomo of the isolated molecule has been cho- 
sen in order to reproduce that calculated by DFT-LDA, 
although variations of ±1 eV around this value do not 
cause any significant change in our results. 

In Fig. 0] we present our calculated I-V curves, the 
occupation of the HOMO and LUMO state as a function 
of bias V, and the DOS for both the weak (<f> = 5°) and 
strong (cf) = 30°) coupling regime. For weak coupling, 
LDA and PSIC give dramatically different I-V charac- 
teristics. In particular, PSIC opens a conductance gap 
in the I-V around zero bias. This is despite that fact 
that the LDA and PSIC DOS look almost identical. In 
both cases E-p pins the bottom of the S-derived empty 
7T* orbital, which is the first state to get involved in the 
transport process. Once bias is applied, such a LUMO 
state gets progressively more populated and follows the 
lead kept at positive bias. The current is then roughly 
proportional to the state occupation, as seen previously 
in the case of the simple model. The key point here is 
that, while in LDA the state charges linearly with bias, 
in PSIC it can follow the upper bias without charging 
significantly. Again the onset of charging will become 
important only when the state has moved upwards in en- 
ergy enough to match the derivative discontinuity. At 
this point the LUMO ir* state starts to conduct (around 
V=l Volt). In addition, for such biases, also the HOMO 
7T state appears in the bias window and contributes to 
the current. 

In the strong coupling limit, the differences between 
LDA and PSIC are much less evident. In this case, both 
the 7T and n* states are very broad (see figure IHb2)) 
providing contributions to the current, even at low bias. 
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FIG. 4: (a) I-V characteristic, (b) occupation as a function 
of energy, and (c) DOS for the system BDT+Au leads. The 
cases of strong (</> = 30°) and weak coupling (<j> = 5°) are 
presented respectively in the right- and left-hand side panels. 
We present only LDA DOS, since PSIC gives almost identical 
results. In panels (c) the vertical lines denotes the position of 
E F . 



We emphasize here that a more rigorous approach 
would be to do, for example, an exact exchange calcu- 
lation [2i|| within NEGF. This is within present com- 
putational capability, but would be costly. Our results 
demonstrate that such a calculation would yield very dif- 
ferent results from LDA or GGA calculations in the weak 
coupling limit . 

In conclusion, we have discussed the main character- 
istics of DFT-based NEGF methods. We have identified 
the lack of derivative discontinuity in continuous density 
functional approximations as a major source of error in 
calculating the I-V characteristic of a molecular junc- 
tion. Our results demonstrate that LDA and GGA are 
not suitable for transport calculations, at least when the 
coupling is weak. We have further proposed a simple 
corrective scheme based on the removal of the atomic 
self-interaction. This has the remarkable property of re- 
introducing, albeit in an approximate way, the derivative 
discontinuity of the potential, while adding moderate ad- 
ditional computational costs. These KS eigenvalues are 
more closely related to the true removal energies, and 
therefore can be employed in a NEGF transport calcu- 
lation. We have implemented such a method in a sim- 
plified tight-bonding scheme and demonstrated that con- 



ductance gaps at low bias can open for molecular junc- 
tions predicted metallic by LDA. 
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